Importing necessary libraries¶
In [1]:
import numpy as np
from scipy.io import wavfile
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d
from scipy.signal import tf2zpk, find_peaks, freqz
from scipy.io.wavfile import write
Please see¶
Sounds are generated and stored in .wav file in the same folder
Please see : Qn7 : /aa/ sound resynthesized - $aaReconstructed.wav$
Qn7 : Pitch variation - $aaReconstructedPitchVariationLow.wav$, $aaReconstructedPitchVariationHigh.wav$
Bonus : /ss/ sound resynthesized - $ssReconstructed.wav$
Question 1¶
Applying pre-emphasis to the signal
We can do this by passing the signal througha filter : $$H(z) = (1 - \alpha z^{-1}) \text{ with $\alpha$ very close to 1, typically = 0.95}$$ $$\implies y[n] = x[n] - \alpha x[n-1]$$
In [2]:
pathToaaSound = r"aa.wav" #path to the "pani" sound of machali
sampleRate, data = wavfile.read(pathToaaSound) #storing the used sample rate and the wav file in variables for analysis
data = data/(32767.0)#Normalising the data
data
Out[2]:
In [3]:
alpha = 0.95#choosing the zero for pre-emphasis
dataPreEmph = np.array([data[n] - alpha*data[n-1] for n in range(1, len(data))])#gives singal from n=1 since we start from 1
dataPreEmph = np.insert(dataPreEmph, 0, data[0])
Question 2¶
Compute and plot the narrowband magnitude spectrum slice using a Hamming window of duration = 30 ms on a segment near the centre of the given audio file.
In [4]:
dur = 0.03 #duration of the window in seconds
windowLength = int(sampleRate*dur) #the window length in samples (each sample is 1/sampleRate long)
totalSignalLength = int(len(dataPreEmph)) #the lenght of the entire signal in terms of samples
signalToAnalyse = dataPreEmph[totalSignalLength//2 - windowLength//2 : totalSignalLength//2 + windowLength//2] #taking the sample of the signal in the center
checkStr = "Window length matches the signal sample Proceed....." if (len(signalToAnalyse)==windowLength) else "Check the signal sample again <Lengths do not match>"
print(checkStr) #checking whether the calculated window lenght and the lenght of windowed signal match
Generating the hamming window¶
I have used the numpy implementation of the hamming window
In [5]:
window = np.hamming(windowLength)
plt.gcf().set_dpi(300)
plt.plot(window)
plt.title("Hamming window")
plt.ylabel("Amplitude")
plt.xlabel("Sample")
plt.show()
Generating the narrawband magnitude spectra¶
Here we multply the hamming window of appropriate length with our signal sample and take its dft. Thus, we can generate the required magnitude plot
In [6]:
dftSize = 10000 #setting the dft size
signalAfterHamming = signalToAnalyse*np.hamming(windowLength) #choosing the hamming window
signalHammingDFT = 20*np.log10(np.abs(np.fft.fft(signalAfterHamming, dftSize))[:dftSize//2]) #computing the magnitudes of the signal afterwindowing and hamming
In [7]:
plt.figure(figsize=(12, 8)) #plotting the magnitude spectra for the windowed signal
plt.gcf().set_dpi(300)
plt.plot(np.arange(0, 4000, 8000/dftSize)[0:] , signalHammingDFT)
plt.xlabel("Frequency in Hertz -->")
plt.ylabel("Magnitude of the DFT in dB -->")
plt.legend(["Magnitude of the DFT"])
Out[7]:
Question 3¶
With the same 30 ms segment of part 2, compute the autocorrelation coefficients required for LPC calculation at various p = 2,4,6,8,10. Use the Levinson-Durbin recursion to compute the LP coefficients from the autocorrelation coefficients. Plot error signal energy (i.e. square of gain) vs p.
For a particular p, we only need to compute the autocorrelation functions from 0 to p-1, since only these coefficients show up in our equations.
To compute the ACF :
$$ r(k) = \sum_{m = -M}^{M} x[m] \times x[m+k]$$$$ \implies r(k) = \sum_{m = 0}^{2M+1} x[m] \times x[m+k] \text{(when signal is indexed from 0)} $$Further, since our signal is windowed i.e finite length the autocorrelation is always 0 for all $ k \geq N \text{(the window length)}$
In [8]:
def computeACF(signal, lag):
#we assume that the signal is windowed
#padding the signal till twice the length since the signal is finite
signal = np.array(signal, dtype = float)
return (np.sum(np.pad(signal, (0, 3*len(signal)), 'constant')*np.pad(signal, (lag, 3*len(signal)-lag), 'constant')))
In [9]:
k = [i for i in range(0, 240)]
x = np.linspace(0, 239, 2000)
acf = [computeACF(signalAfterHamming, i) for i in range(0, 240)]
plt.figure(figsize = (20, 8))
plt.gcf().set_dpi(300)
plt.plot(k, acf,'.')
f = interp1d(k, acf, kind='cubic')
plt.plot(x, f(x), alpha = 0.3, linewidth = 0.9, color = 'g')
plt.vlines(k, 0, acf, color='r', linestyles='--', alpha = 0.3)
plt.axhline(0, color = 'k', alpha = 0.4)
plt.legend(["ACF at different lag values", "Interpolated ACF plot"])
plt.title("AutoCorrelation function at different lag values", fontsize= 15)
plt.xlabel("Lag (k)", fontsize= 15)
plt.ylabel("ACF Value", fontsize= 15)
Out[9]:
In [10]:
plt.figure(figsize=(12, 16))
plt.suptitle("Signal of the sound at different positions")
plt.subplot(311)
lag = 0
plt.plot(np.pad(signalAfterHamming,(lag, 3*len(signalAfterHamming)-lag), 'constant'))
plt.xlabel("Samples[n]")
plt.ylabel("Amplitude")
plt.title("Original Signal")
plt.subplot(312)
lag = 200
plt.plot(np.pad(signalAfterHamming,(lag, 3*len(signalAfterHamming)-lag), 'constant'))
plt.xlabel("Samples[n]")
plt.ylabel("Amplitude")
plt.title(f"Shifted Signal by {lag}")
plt.subplot(313)
lag = 200
plt.plot(np.pad(signalAfterHamming,(0, 3*len(signalAfterHamming)), 'constant')*np.pad(signalAfterHamming,(lag, 3*len(signalAfterHamming)-lag), 'constant'))
plt.xlabel("Samples[n]")
plt.ylabel("Amplitude")
plt.title(f"Signal product actually obtained at {lag}")
Out[10]:
The Levinson Durban Recursion¶
We use the LD recursion as mentioned in class.
Emin is computed using the autocorrelation values
$$Emin =\left( r[0] - \sum_{k=1}^{P} (r[k] \times a_k) \right) \text{\;where P is the order of LD recursion used}$$$$Error = Emin = \sum_{n} \left ( s[n] - \sum_{k=1}^{P} (a_k \times s[n-k]) \right)^2$$In [11]:
rp = np.array([computeACF(signalAfterHamming, lag = 1)])
pp = np.array([computeACF(signalAfterHamming, lag = 1)])
alpha = np.array([rp[0]/computeACF(signalAfterHamming, 0)])
beta = np.array([rp[0]/computeACF(signalAfterHamming, 0)])
coef = [alpha]
order = 40
for p in range(1, order):
k = (computeACF(signalAfterHamming, p+1) - np.sum(pp*alpha))/(computeACF(signalAfterHamming, 0) - np.sum(rp*alpha))
alpha = np.append(alpha, 0) - (k*np.append(beta, -1))
beta = alpha[::-1]
rp = np.append(rp, computeACF(signalAfterHamming, lag = p+1))
pp = rp[::-1]
coef.append(alpha)
In [12]:
coef
Out[12]:
In [13]:
#checking whether the coefficient are correct
#should preserve the acf!
for p in range(1, len(coef)+1):
computeACF(signalAfterHamming, 1) == coef[0][0]*computeACF(signalAfterHamming, 0)
Rind = np.array([[ np.abs(i-k) for k in range(1, p+1)] for i in range(1, p+1)])
R = np.array([[computeACF(signalAfterHamming, np.abs(i-k)) for k in range(1, p+1)] for i in range(1, p+1)])
acfArr = np.array([computeACF(signalAfterHamming, i) for i in range(1, p+1)])
coefArr = coef[p-1]
if np.sum((np.round(np.matmul(R, coefArr) - acfArr, 7))**2)==0:
print(f"Order {p} coefficients verified and correct...proceed")
else:
print("Error <Coefficients don't match>")
break
In [14]:
#computing the minima energy directly from formula
Emin = []
for p in range(1,len(coef)+1):
acfArr = np.array([computeACF(signalAfterHamming, i) for i in range(1, p+1)])
E = computeACF(signalAfterHamming, 0) - np.sum(acfArr*coef[p-1])
Emin.append(E)
In [15]:
error = []
for p in range(1, len(coef)+1):
signalShifted = np.pad(np.array(signalAfterHamming, dtype = float),(p, 3*len(signalAfterHamming)-p), 'constant')
signalToEstimate = np.pad(np.array(signalAfterHamming, dtype = float),(0, 3*len(signalAfterHamming)), 'constant')
signalEstimate = np.zeros_like(signalShifted)
for i in range(0, len(signalAfterHamming) + p-1):
signalEstimate[i] = np.sum(signalShifted[i:(i+p)]*coef[p-1][::-1])
error.append(np.sum((signalToEstimate - signalEstimate)**2))
Error plots¶
In [16]:
plt.figure(figsize=(16, 8))
plt.gcf().set_dpi(300)
plt.suptitle("Error computed for different orders of LPA for different methods of error computation")
plt.subplot(121)
plt.plot([i for i in range(1, order + 1)], error, 'ro')
plt.plot([i for i in range(1, order + 1)], error, 'b', alpha = 0.3, linewidth = 1)
plt.vlines([i for i in range(1, order + 1)] , 0, error, color='g', linestyles='--', alpha = 0.3)
plt.ylim(bottom=0)
plt.title("Calculated MSE error/loss from the computed coefficients")
plt.xlabel("Order of LPA (p)")
plt.ylabel("Calculated error")
plt.legend(["Error computed"])
plt.subplot(122)
plt.plot([i for i in range(1, order + 1)], Emin, 'ro')
plt.plot([i for i in range(1, order + 1)], Emin, 'b', alpha = 0.3, linewidth = 1)
plt.vlines([i for i in range(1, order + 1)] , 0, Emin, color='g', linestyles='--', alpha = 0.3)
plt.ylim(bottom=0)
plt.title("Calculated MSE error/loss from direct loss equation")
plt.xlabel("Order of LPA (p)")
plt.ylabel("Cmputed error from formula")
plt.legend(["Error computed from formula"])
Out[16]:
Inference¶
- We see that the computed loss reduce as p increases and finally settle at minima around p=40
- This is expected since higher order lpa approximated the system better than lower order lpa
Question 4¶
Show the pole-zero plots of the estimated all-pole filter for p=6,10; comment.
Order of lpa p=6¶
In [17]:
z, p, k =tf2zpk(np.sqrt(Emin[5]), np.insert((-1)*coef[5] , 0, 1))
In [18]:
#Zeros
z
Out[18]:
In [19]:
#poles
p
Out[19]:
In [20]:
angle = np.linspace( 0 , 2 * np.pi , 150 )
radius = 1
x = radius * np.cos( angle )
y = radius * np.sin( angle )
plt.figure(figsize=(10,10))
plt.plot(x, y, 'r', alpha = 0.2)
plt.plot([i.real for i in p], [j.imag for j in p], 'x', label= "poles")
plt.xlim([-1.1, 1.1])
plt.ylim([-1.1, 1.1])
plt.title("Pole zero plot for LPA of order 6")
plt.xlabel("Real axis")
plt.ylabel("Imaginary axis")
plt.legend()
Out[20]:
Order of lpa p=10¶
In [21]:
z, p, k =tf2zpk(np.sqrt(Emin[9]),np.insert((-1)*coef[9] , 0, 1))
In [22]:
#Zeros
z
Out[22]:
In [23]:
#poles
p
Out[23]:
In [24]:
angle = np.linspace( 0 , 2 * np.pi , 150 )
radius = 1
x = radius * np.cos( angle )
y = radius * np.sin( angle )
plt.figure(figsize=(10,10))
plt.plot(x, y, 'r', alpha = 0.2)
plt.plot([i.real for i in p], [j.imag for j in p], 'x', label = "Poles")
plt.xlim([-1.1, 1.1])
plt.ylim([-1.1, 1.1])
plt.title("Pole zero plot for LPA of order 10")
plt.xlabel("Real axis")
plt.ylabel("Imaginary axis")
plt.legend()
Out[24]:
Comments:¶
- We see that most poles lie very close to the unit circle in both the cases(for p = 6 a few poles are a little inside the unit circle)
- We see 2 poles lying on the real axis in the case of p=10 which is unexpected!
Question 5¶
- Compute the gain
and - plot the LPC spectrum magnitude (i.e. the dB magnitude frequency response of the estimated all-pole filter) for each order "p".
- Comment on the characteristics of the spectral envelope estimates.
- Comment on their shapes with reference to the short-time magnitude spectrum computed in part 2.
Gain¶
Gain is obtained as follows $G = \sqrt{E_{min}}$
In [25]:
#can we obtained by sqrt(Emin)
for i in range(0, order):
print(f"Gain for order {i+1} lpa is {np.sqrt(Emin[i])}")
LPC spectra magnitude plots¶
For order p = 2¶
In [26]:
p = 2
signalToEstimate = np.pad(np.array(signalToAnalyse, dtype = float),(p, 3*len(signalToAnalyse)-p), 'constant')
signalEstimate = np.zeros_like(signalToEstimate)
for i in range(0, len(signalToAnalyse) + p-1):
signalEstimate[i] = np.sum(signalToEstimate[i:(i+p)]*coef[p-1])
w, h = freqz(np.sqrt(Emin[p-1]), np.insert((-1*coef[p-1]), 0, 1), None, 1)
plt.figure(figsize=(12, 8)) #plotting the magnitude spectra for the windowed signal
plt.gcf().set_dpi(300)
#plt.subplot(121)
plt.plot(np.arange(0, 8000/2, 8000/dftSize)[0:], 20*np.log10(np.abs(np.fft.fft(signalAfterHamming, dftSize))[:dftSize//2]))
plt.plot(w[:int(len(w)/2)]*8000/(2*3.14), 20*np.log10(np.abs(h[:int(len(h)/2)])))
plt.xlabel("Frequency in Hertz -->")
plt.ylabel("Magnitude of the DFT in dB -->")
plt.legend(["Magnitude of the DFT", f"Magnitude of the LPC of order {p}"])
# plt.subplot(122)
# f = interp1d(find_peaks(20*np.log(np.abs(np.fft.fft(signalToEstimate, dftSize))[:dftSize//2]))[0], 20*np.log(np.abs(np.fft.fft(signalToEstimate, dftSize))[:dftSize//2])[find_peaks(20*np.log(np.abs(np.fft.fft(signalToEstimate, dftSize))[:dftSize//2]))[0]], kind='cubic')
# xNew = np.array([i for i in range(21, 4000)])
# plt.plot(xNew, f(xNew), alpha = 0.3, linewidth = 0.9, color = 'g')
# plt.xlabel("Frequency in Hertz -->")
# plt.ylabel("Magnitude of the DFT in dB -->")
# plt.legend(["Interpolated Magnitude of the DFT"])
Out[26]:
For order p = 4¶
In [27]:
p = 4
signalToEstimate = np.pad(np.array(signalToAnalyse, dtype = float),(p, 3*len(signalToAnalyse)-p), 'constant')
signalEstimate = np.zeros_like(signalToEstimate)
for i in range(0, len(signalToAnalyse) + p-1):
signalEstimate[i] = np.sum(signalToEstimate[i:(i+p)]*coef[p-1])
w, h = freqz(np.sqrt(Emin[p-1]), np.insert((-1*coef[p-1]), 0, 1), None, 1)
plt.figure(figsize=(12, 8)) #plotting the magnitude spectra for the windowed signal
plt.gcf().set_dpi(300)
#plt.subplot(121)
plt.plot(np.arange(0, 8000/2, 8000/dftSize)[0:], 20*np.log10(np.abs(np.fft.fft(signalAfterHamming, dftSize))[:dftSize//2]))
plt.plot(w[:int(len(w)/2)]*8000/(2*3.14), 20*np.log10(np.abs(h[:int(len(h)/2)])))
plt.xlabel("Frequency in Hertz -->")
plt.ylabel("Magnitude of the DFT in dB -->")
plt.legend(["Magnitude of the DFT", f"Magnitude of the LPC of order {p}"])
# plt.subplot(122)
# f = interp1d(find_peaks(20*np.log(np.abs(np.fft.fft(signalToEstimate, dftSize))[:dftSize//2]))[0], 20*np.log(np.abs(np.fft.fft(signalToEstimate, dftSize))[:dftSize//2])[find_peaks(20*np.log(np.abs(np.fft.fft(signalToEstimate, dftSize))[:dftSize//2]))[0]], kind='cubic')
# xNew = np.array([i for i in range(21, 4000)])
# plt.plot(xNew, f(xNew), alpha = 0.3, linewidth = 0.9, color = 'g')
# plt.xlabel("Frequency in Hertz -->")
# plt.ylabel("Magnitude of the DFT in dB -->")
# plt.legend(["Interpolated Magnitude of the DFT"])
Out[27]:
For order p = 6¶
In [28]:
p = 6
signalToEstimate = np.pad(np.array(signalToAnalyse, dtype = float),(p, 3*len(signalToAnalyse)-p), 'constant')
signalEstimate = np.zeros_like(signalToEstimate)
for i in range(0, len(signalToAnalyse) + p-1):
signalEstimate[i] = np.sum(signalToEstimate[i:(i+p)]*coef[p-1])
w, h = freqz(np.sqrt(Emin[p-1]), np.insert((-1*coef[p-1]), 0, 1), None, 1)
plt.figure(figsize=(12, 8)) #plotting the magnitude spectra for the windowed signal
plt.gcf().set_dpi(300)
#plt.subplot(121)
plt.plot(np.arange(0, 8000/2, 8000/dftSize)[0:], 20*np.log10(np.abs(np.fft.fft(signalAfterHamming, dftSize))[:dftSize//2]))
plt.plot(w[:int(len(w)/2)]*8000/(2*3.14), 20*np.log10(np.abs(h[:int(len(h)/2)])))
plt.xlabel("Frequency in Hertz -->")
plt.ylabel("Magnitude of the DFT in dB -->")
plt.legend(["Magnitude of the DFT", f"Magnitude of the LPC of order {p}"])
# plt.subplot(122)
# f = interp1d(find_peaks(20*np.log(np.abs(np.fft.fft(signalToEstimate, dftSize))[:dftSize//2]))[0], 20*np.log(np.abs(np.fft.fft(signalToEstimate, dftSize))[:dftSize//2])[find_peaks(20*np.log(np.abs(np.fft.fft(signalToEstimate, dftSize))[:dftSize//2]))[0]], kind='cubic')
# xNew = np.array([i for i in range(21, 4000)])
# plt.plot(xNew, f(xNew), alpha = 0.3, linewidth = 0.9, color = 'g')
# plt.xlabel("Frequency in Hertz -->")
# plt.ylabel("Magnitude of the DFT in dB -->")
# plt.legend(["Interpolated Magnitude of the DFT"])
Out[28]:
For order p = 8¶
In [29]:
p = 8
signalToEstimate = np.pad(np.array(signalToAnalyse, dtype = float),(p, 3*len(signalToAnalyse)-p), 'constant')
signalEstimate = np.zeros_like(signalToEstimate)
for i in range(0, len(signalToAnalyse) + p-1):
signalEstimate[i] = np.sum(signalToEstimate[i:(i+p)]*coef[p-1])
w, h = freqz(np.sqrt(Emin[p-1]), np.insert((-1*coef[p-1]), 0, 1), None, 1)
plt.figure(figsize=(12, 8)) #plotting the magnitude spectra for the windowed signal
plt.gcf().set_dpi(300)
#plt.subplot(121)
plt.plot(np.arange(0, 8000/2, 8000/dftSize)[0:], 20*np.log10(np.abs(np.fft.fft(signalAfterHamming, dftSize))[:dftSize//2]))
plt.plot(w[:int(len(w)/2)]*8000/(2*3.14), 20*np.log10(np.abs(h[:int(len(h)/2)])))
plt.xlabel("Frequency in Hertz -->")
plt.ylabel("Magnitude of the DFT in dB -->")
plt.legend(["Magnitude of the DFT", f"Magnitude of the LPC of order {p}"])
# plt.subplot(122)
# f = interp1d(find_peaks(20*np.log(np.abs(np.fft.fft(signalToEstimate, dftSize))[:dftSize//2]))[0], 20*np.log(np.abs(np.fft.fft(signalToEstimate, dftSize))[:dftSize//2])[find_peaks(20*np.log(np.abs(np.fft.fft(signalToEstimate, dftSize))[:dftSize//2]))[0]], kind='cubic')
# xNew = np.array([i for i in range(21, 4000)])
# plt.plot(xNew, f(xNew), alpha = 0.3, linewidth = 0.9, color = 'g')
# plt.xlabel("Frequency in Hertz -->")
# plt.ylabel("Magnitude of the DFT in dB -->")
# plt.legend(["Interpolated Magnitude of the DFT"])
Out[29]:
For order p = 10¶
In [30]:
p = 10
signalToEstimate = np.pad(np.array(signalToAnalyse, dtype = float),(p, 3*len(signalToAnalyse)-p), 'constant')
signalEstimate = np.zeros_like(signalToEstimate)
for i in range(0, len(signalToAnalyse) + p-1):
signalEstimate[i] = np.sum(signalToEstimate[i:(i+p)]*coef[p-1])
w, h = freqz(np.sqrt(Emin[p-1]), np.insert((-1*coef[p-1]), 0, 1), None, 1)
plt.figure(figsize=(12, 8)) #plotting the magnitude spectra for the windowed signal
plt.gcf().set_dpi(300)
#plt.subplot(121)
plt.plot(np.arange(0, 8000/2, 8000/dftSize)[0:], 20*np.log10(np.abs(np.fft.fft(signalAfterHamming, dftSize))[:dftSize//2]))
plt.plot(w[:int(len(w)/2)]*8000/(2*3.14), 20*np.log10(np.abs(h[:int(len(h)/2)])))
plt.xlabel("Frequency in Hertz -->")
plt.ylabel("Magnitude of the DFT in dB -->")
plt.legend(["Magnitude of the DFT", f"Magnitude of the LPC of order {p}"])
# plt.subplot(122)
# f = interp1d(find_peaks(20*np.log(np.abs(np.fft.fft(signalToEstimate, dftSize))[:dftSize//2]))[0], 20*np.log(np.abs(np.fft.fft(signalToEstimate, dftSize))[:dftSize//2])[find_peaks(20*np.log(np.abs(np.fft.fft(signalToEstimate, dftSize))[:dftSize//2]))[0]], kind='cubic')
# xNew = np.array([i for i in range(21, 4000)])
# plt.plot(xNew, f(xNew), alpha = 0.3, linewidth = 0.9, color = 'g')
# plt.xlabel("Frequency in Hertz -->")
# plt.ylabel("Magnitude of the DFT in dB -->")
# plt.legend(["Interpolated Magnitude of the DFT"])
Out[30]:
We clearly see that as the order increases of the lpa analysis we get our estimated envelope more and more closer to the actual waveform spectra
This means that we are able to better identfy the formants (since hidden envelope is estimated) as the order of the lpa increases
Question 6¶
Based on the 10th-order LP coefficients, carry out the inverse filtering of the /a/ vowel segment to obtain the residual error signal. Can you measure the pitch period of the voiced sound from the residual waveform? Use the acf to detect the pitch. Compare the acf plots of the original speech and residual signals.
In [31]:
p = 10
signalToEstimate = np.pad(np.array(signalAfterHamming, dtype = float),(p, 3*len(signalAfterHamming)-p), 'constant')
signalEstimate = np.zeros_like(signalToEstimate)
for i in range(0, len(signalAfterHamming) + p-1):
signalEstimate[i] = np.sum(signalToEstimate[i:(i+p)]*coef[p-1])
In [32]:
plt.figure(figsize = (16, 8))
plt.gcf().set_dpi(300)
plt.suptitle(f"Plots of the signal, its estimate and the error obtained in time domain for order {p}")
plt.subplot(121)
plt.plot(signalToEstimate[p:240+p], label = "Signal(Ground truth)")
plt.plot(signalEstimate[0:240], label = "Signal estimated")
plt.title("Signal along with its estimate")
plt.xlabel("Time(sec)")
plt.ylabel("Amplitude")
plt.legend()
plt.subplot(122)
plt.plot((signalEstimate -signalToEstimate)[0:240], label = "Error")
plt.title("Error in the estimate")
plt.xlabel("Time(sec)")
plt.ylabel("Amplitude")
plt.legend()
plt.legend()
Out[32]:
In [33]:
k = [i for i in range(0, 300)]
x = np.linspace(0, 299, 2000)
acf = [computeACF(signalEstimate -signalToEstimate, i) for i in range(0, 300)]
plt.figure(figsize = (20, 8))
plt.gcf().set_dpi(300)
plt.plot(k, acf,'.')
f = interp1d(k, acf, kind='cubic')
plt.plot(x, f(x), alpha = 0.3, linewidth = 0.9, color = 'g')
plt.vlines(k, 0, acf, color='r', linestyles='--', alpha = 0.3)
plt.axhline(0, color = 'k', alpha = 0.4)
plt.legend(["ACF at different lag values", "Interpolated ACF plot"])
plt.title("AutoCorrelation function at different lag values of the error signal", fontsize= 15)
plt.xlabel("Lag (k)", fontsize= 15)
plt.ylabel("ACF Value", fontsize= 15)
Out[33]:
Referring the Question 3 we have the original ACF plot of the signal.
Comparing the two, we have :
- Both of them peak at nearly 59 samples (i.e 59 $\times \frac{1}{8000} = 7.375$ micro-seconds).
- The above behaviour is expected since, both signals have the underlying pitch(glottal vabrations)
- Thus, the estimated pitch is $\frac{1}{7.375ms} = 135.925Hz$
- The nature of the two ACF graphs is different (i.e.) their shape etc and this is obvious since the underlying signal is different but still point (1) holds
Question 7¶
LP re-synthesis: We analysed a natural speech sound /a/ above.
Using a suitable set of parameter estimates as obtained there, we wish to reconstruct the sound. That is, use the best estimated LP filter with an ideal impulse train of the estimated pitch period as source excitation.
Carry out de-emphasis on the output waveform. Set the duration of the synthesized sound to be 300 ms at 8 kHz sampling frequency and view the waveform as well as listen to your created sound.
Comment on the similarity with the original sound. Try out voice modification using this analysis-synthesis method (e.g. change the voice pitch).
We will use the order 10 filter.
Say u[n] is our impulse train. Clearly, u[n] = 1 when n = $k \times T$ else u[n] = 0 where T is the pitch of the signal i.e. 59 samples as calculated in the previous question.
To calculate the output(estimated signal), we need to set up the difference equation as follows:
We know our filter looks like : $$A(z) = \frac{G}{1 - \sum_{k = 1}^{k=P}a_k z^{-k}}$$ $$\implies y[n] - \sum_{k = 1}^{k=P}a_k y[n-k] = Gx[n]$$
$$\implies y[n] = Gx[n] + \sum_{k = 1}^{k=P}a_k y[n-k]$$We assume that the signal y[n] is causal and x[n] is our impulse i.e. impulse train
In [34]:
signalToEstimate = np.pad(np.array(signalAfterHamming, dtype = float),(0, 10*len(signalAfterHamming)), 'constant')
signalReconstructed = np.zeros_like(signalToEstimate)
inputImpulseTrain = np.zeros_like(signalToEstimate)
inputImpulseTrain[[i for i in range(0, signalToEstimate.shape[0], 59)]] = 1
#inputImpulseTrain[0]=1
plt.figure(figsize = (20, 8))
plt.gcf().set_dpi(300)
plt.plot(inputImpulseTrain, 'o')
plt.vlines(np.arange(0, signalToEstimate.shape[0], 59), 0, 1, colors = 'r', linestyles='--')
plt.xlabel("Times(samples)")
plt.ylabel("Amplitude")
plt.title("Impulse train generated")
Out[34]:
In [35]:
P = 10#choosing the order of analysis
G = np.sqrt(Emin[P-1])
inputImpulseTrain = np.pad(np.array(inputImpulseTrain, dtype = float),(P, P), 'constant', constant_values=0)
#signal y[n] is assumed casual
for i in range(P, signalToEstimate.shape[0]):
signalReconstructed[i] = G*inputImpulseTrain[i] + np.sum(signalReconstructed[i-P:i]*coef[P-1][::-1])
We now need to de-emphasis, we can do this by passing the signal through a filter : $$H(z) = \frac{1}{(1 - \alpha z^{-1})} \text{ with $\alpha$ very close to 1, typically = 0.95}$$ $$\implies s_{de-emph}[n] = s_{emph}[n] + \alpha s_{de-emph}[n-1]$$
In [36]:
alpha = 0.95#setting the same value as pre-emphasis
signalReconstructed = signalReconstructed[P:]
signalReconstructedDeEmph = np.pad(np.zeros_like(signalReconstructed), (1, 0), 'constant')
for j in range(1, signalReconstructedDeEmph.shape[0]-1):
signalReconstructedDeEmph[j] = signalReconstructed[j] + alpha*signalReconstructedDeEmph[j-1]
In [37]:
plt.figure(figsize=(16,8))
plt.gcf().set_dpi(300)
plt.suptitle("Signal estimated along with ground truth")
plt.subplot(121)
plt.plot([i*(1.0/8000) for i in range(len(signalReconstructed))][0:240],signalReconstructedDeEmph[1:][0:240])
plt.xlabel("Time(seconds)")
plt.ylabel("Amplitude")
plt.title("Estimate of the signal(after resynthesis and de-emphasis)")
plt.subplot(122)
plt.plot([i*(1.0/8000) for i in range(len(signalReconstructed))][0:240],signalToAnalyse)
plt.xlabel("Time(seconds)")
plt.ylabel("Amplitude")
plt.title("Ground truth(true signal")
Out[37]:
In [38]:
scaled = np.int16(signalReconstructedDeEmph/np.max(np.abs(signalReconstructedDeEmph)) * 32767)#scaling the output to 16bit integers for writing to wav file
write('aaReconstructed.wav', 8000, scaled)
We see the sound is very similar to the original sound
Further the waveforms are similar looking -notice their amplitudes etc. may not match exaclty
Changing pitch of the sound¶
Notice that the beauty of the lpa analysis is that we are finding the envelope of the sound that is the glottal vibrations are not taken into account. This means that all information about the glottal vibrations that is the pitch is encoded in the e[n] signal that is fed into our model. Here e[n] is the ideal impulse train where the distance(in samples) between two impulses encodes the information about the pitch. Thus, modifying the impulse train basically gives us a rough estimate of the sound at a different pitch
Changing pitch to higher pitch¶
We would need to reduce the T0 between each impulse thus F0 increases and so does Pitch
In [39]:
signalToEstimate = np.pad(np.array(signalAfterHamming, dtype = float),(0, 10*len(signalAfterHamming)), 'constant')
signalReconstructed = np.zeros_like(signalToEstimate)
inputImpulseTrain = np.zeros_like(signalToEstimate)
inputImpulseTrain[[i for i in range(0, signalToEstimate.shape[0], 20)]] = 1
#inputImpulseTrain[0]=1
plt.figure(figsize = (20, 8))
plt.gcf().set_dpi(300)
plt.plot(inputImpulseTrain, 'o')
plt.vlines(np.arange(0, signalToEstimate.shape[0], 20), 0, 1, colors = 'r', linestyles='--')
plt.xlabel("Times(samples)")
plt.ylabel("Amplitude")
plt.title("Impulse train generated")
Out[39]:
In [40]:
P = 10#choosing the order of analysis
G = np.sqrt(Emin[P-1])
inputImpulseTrain = np.pad(np.array(inputImpulseTrain, dtype = float),(P, P), 'constant', constant_values=0)
#signal y[n] is assumed casual
for i in range(P, signalToEstimate.shape[0]):
signalReconstructed[i] = G*inputImpulseTrain[i] + np.sum(signalReconstructed[i-P:i]*coef[P-1][::-1])
In [41]:
alpha = 0.95#setting the same value as pre-emphasis
signalReconstructed = signalReconstructed[P:]
signalReconstructedDeEmph = np.pad(np.zeros_like(signalReconstructed), (1, 0), 'constant')
for j in range(1, signalReconstructedDeEmph.shape[0]-1):
signalReconstructedDeEmph[j] = signalReconstructed[j] + alpha*signalReconstructedDeEmph[j-1]
In [42]:
plt.figure(figsize=(16,8))
plt.gcf().set_dpi(300)
plt.suptitle("Signal estimated along with ground truth")
plt.subplot(121)
plt.plot([i*(1.0/8000) for i in range(len(signalReconstructed))][0:240],signalReconstructedDeEmph[1:][0:240])
plt.xlabel("Time(seconds)")
plt.ylabel("Amplitude")
plt.title("Estimate of the signal(after resynthesis and de-emphasis)")
plt.subplot(122)
plt.plot([i*(1.0/8000) for i in range(len(signalReconstructed))][0:240],signalToAnalyse)
plt.xlabel("Time(seconds)")
plt.ylabel("Amplitude")
plt.title("Estimate of the signal(after resynthesis)")
Out[42]:
In [43]:
scaled = np.int16(signalReconstructedDeEmph/np.max(np.abs(signalReconstructedDeEmph)) * 32767)#scaling the output to 16bit integers for writing to wav file
write('aaReconstructedPitchVariationHigh.wav', 8000, scaled)
Changing pitch to lower pitch¶
We would need to increase the T0 between each impulse thus F0 decreases and so does Pitch
In [44]:
signalToEstimate = np.pad(np.array(signalAfterHamming, dtype = float),(0, 10*len(signalAfterHamming)), 'constant')
signalReconstructed = np.zeros_like(signalToEstimate)
inputImpulseTrain = np.zeros_like(signalToEstimate)
inputImpulseTrain[[i for i in range(0, signalToEstimate.shape[0], 100)]] = 1
#inputImpulseTrain[0]=1
plt.figure(figsize = (20, 8))
plt.gcf().set_dpi(300)
plt.plot(inputImpulseTrain, 'o')
plt.vlines(np.arange(0, signalToEstimate.shape[0], 100), 0, 1, colors = 'r', linestyles='--')
plt.xlabel("Times(samples)")
plt.ylabel("Amplitude")
plt.title("Impulse train generated")
Out[44]:
In [45]:
P = 10#choosing the order of analysis
G = np.sqrt(Emin[P-1])
inputImpulseTrain = np.pad(np.array(inputImpulseTrain, dtype = float),(P, P), 'constant', constant_values=0)
#signal y[n] is assumed casual
for i in range(P, signalToEstimate.shape[0]):
signalReconstructed[i] = G*inputImpulseTrain[i] + np.sum(signalReconstructed[i-P:i]*coef[P-1][::-1])
In [46]:
alpha = 0.95#setting the same value as pre-emphasis
signalReconstructed = signalReconstructed[P:]
signalReconstructedDeEmph = np.pad(np.zeros_like(signalReconstructed), (1, 0), 'constant')
for j in range(1, signalReconstructedDeEmph.shape[0]-1):
signalReconstructedDeEmph[j] = signalReconstructed[j] + alpha*signalReconstructedDeEmph[j-1]
In [47]:
plt.figure(figsize=(16,8))
plt.gcf().set_dpi(300)
plt.suptitle("Signal estimated along with ground truth(shown for few seconds for better clarity)")
plt.subplot(121)
plt.plot([i*(1.0/8000) for i in range(len(signalReconstructed))][0:240],signalReconstructedDeEmph[1:][0:240])
plt.xlabel("Time(seconds)")
plt.ylabel("Amplitude")
plt.title("Estimate of the signal(after resynthesis and de-emphasis)")
plt.subplot(122)
plt.plot([i*(1.0/8000) for i in range(len(signalReconstructed))][0:240],signalToAnalyse)
plt.xlabel("Time(seconds)")
plt.ylabel("Amplitude")
plt.title("Estimate of the signal(after resynthesis)")
Out[47]:
In [48]:
scaled = np.int16(signalReconstructedDeEmph/np.max(np.abs(signalReconstructedDeEmph)) * 32767)#scaling the output to 16bit integers for writing to wav file
write('aaReconstructedPitchVariationLow.wav', 8000, scaled)
[Bonus]¶
Perform LP analysis on the provided /s/ sound (sampled at 16 kHz).
Same steps are repeated as above
In [49]:
pathToaaSound = r"ss.wav" #path to the "pani" sound of machali
sampleRate, data = wavfile.read(pathToaaSound) #storing the used sample rate and the wav file in variables for analysis
data = data/(32767.0)
data
Out[49]:
In [50]:
dur = 0.03 #duration of the window in seconds
windowLength = int(sampleRate*dur) #the window length in samples (each sample is 1/sampleRate long)
totalSignalLength = int(len(data)) #the lenght of the entire signal in terms of samples
signalToAnalyse = data[totalSignalLength//2 - windowLength//2 : totalSignalLength//2 + windowLength//2] #taking the sample of the signal in the center
checkStr = "Window length matches the signal sample Proceed....." if (len(signalToAnalyse)==windowLength) else "Check the signal sample again <Lengths do not match>"
print(checkStr) #checking whether the calculated window lenght and the lenght of windowed signal match
In [51]:
dftSize = 10000 #setting the dft size
signalAfterHamming = signalToAnalyse*np.hamming(windowLength) #choosing the hamming window
signalHammingDFT = 20*np.log10(np.abs(np.fft.fft(signalAfterHamming, dftSize))[:dftSize//2]) #computing the magnitudes of the signal afterwindowing and hamming
In [52]:
plt.figure(figsize=(12, 8)) #plotting the magnitude spectra for the windowed signal
plt.gcf().set_dpi(300)
plt.plot(np.arange(0, 4000, 8000/dftSize)[0:] , signalHammingDFT)
plt.xlabel("Frequency in Hertz -->")
plt.ylabel("Magnitude of the DFT in dB -->")
plt.legend(["Magnitude of the DFT"])
Out[52]:
In [53]:
k = [i for i in range(0, 240)]
x = np.linspace(0, 239, 2000)
acf = [computeACF(signalAfterHamming, i) for i in range(0, 240)]
plt.figure(figsize = (20, 8))
plt.gcf().set_dpi(300)
plt.plot(k, acf,'.')
f = interp1d(k, acf, kind='cubic')
plt.plot(x, f(x), alpha = 0.3, linewidth = 0.9, color = 'g')
plt.vlines(k, 0, acf, color='r', linestyles='--', alpha = 0.3)
plt.axhline(0, color = 'k', alpha = 0.4)
plt.legend(["ACF at different lag values", "Interpolated ACF plot"])
plt.title("AutoCorrelation function at different lag values", fontsize= 15)
plt.xlabel("Lag (k)", fontsize= 15)
plt.ylabel("ACF Value", fontsize= 15)
Out[53]:
In [54]:
rp = np.array([computeACF(signalAfterHamming, lag = 1)])
pp = np.array([computeACF(signalAfterHamming, lag = 1)])
alpha = np.array([rp[0]/computeACF(signalAfterHamming, 0)])
beta = np.array([rp[0]/computeACF(signalAfterHamming, 0)])
coef = [alpha]
order = 40
for p in range(1, order):
k = (computeACF(signalAfterHamming, p+1) - np.sum(pp*alpha))/(computeACF(signalAfterHamming, 0) - np.sum(rp*alpha))
alpha = np.append(alpha, 0) - (k*np.append(beta, -1))
beta = alpha[::-1]
rp = np.append(rp, computeACF(signalAfterHamming, lag = p+1))
pp = rp[::-1]
coef.append(alpha)
In [55]:
coef
Out[55]:
In [56]:
#checking whether the coefficient are correct
#should preserve the acf!
for p in range(1, len(coef)+1):
computeACF(signalAfterHamming, 1) == coef[0][0]*computeACF(signalAfterHamming, 0)
Rind = np.array([[ np.abs(i-k) for k in range(1, p+1)] for i in range(1, p+1)])
R = np.array([[computeACF(signalAfterHamming, np.abs(i-k)) for k in range(1, p+1)] for i in range(1, p+1)])
acfArr = np.array([computeACF(signalAfterHamming, i) for i in range(1, p+1)])
coefArr = coef[p-1]
if np.sum((np.round(np.matmul(R, coefArr) - acfArr, 7))**2)==0:
print(f"Order {p} coefficients verified and correct...proceed")
else:
print("Error <Coefficients don't match>")
break
In [57]:
#computing the minima energy directly from formula
Emin = []
for p in range(1,len(coef)+1):
acfArr = np.array([computeACF(signalAfterHamming, i) for i in range(1, p+1)])
E = computeACF(signalAfterHamming, 0) - np.sum(acfArr*coef[p-1])
Emin.append(E)
In [58]:
error = []
for p in range(1, len(coef)+1):
signalShifted = np.pad(np.array(signalAfterHamming, dtype = float),(p, 3*len(signalAfterHamming)-p), 'constant')
signalToEstimate = np.pad(np.array(signalAfterHamming, dtype = float),(0, 3*len(signalAfterHamming)), 'constant')
signalEstimate = np.zeros_like(signalShifted)
for i in range(0, len(signalAfterHamming) + p-1):
signalEstimate[i] = np.sum(signalShifted[i:(i+p)]*coef[p-1][::-1])
error.append(np.sum((signalToEstimate - signalEstimate)**2))
In [59]:
plt.figure(figsize=(16, 8))
plt.gcf().set_dpi(300)
plt.suptitle("Error computed for different orders of LPA for different methods of error computation")
plt.subplot(121)
plt.plot([i for i in range(1, order + 1)], error, 'ro')
plt.plot([i for i in range(1, order + 1)], error, 'b', alpha = 0.3, linewidth = 1)
plt.vlines([i for i in range(1, order + 1)] , 0, error, color='g', linestyles='--', alpha = 0.3)
plt.ylim(bottom=0)
plt.title("Calculated MSE error/loss from the computed coefficients")
plt.xlabel("Order of LPA (p)")
plt.ylabel("Calculated error")
plt.legend(["Error computed"])
plt.subplot(122)
plt.plot([i for i in range(1, order + 1)], Emin, 'ro')
plt.plot([i for i in range(1, order + 1)], Emin, 'b', alpha = 0.3, linewidth = 1)
plt.vlines([i for i in range(1, order + 1)] , 0, Emin, color='g', linestyles='--', alpha = 0.3)
plt.ylim(bottom=0)
plt.title("Calculated MSE error/loss from direct loss equation")
plt.xlabel("Order of LPA (p)")
plt.ylabel("Cmputed error from formula")
plt.legend(["Error computed from formula"])
Out[59]:
Order of lpa p=6¶
In [60]:
z, p, k =tf2zpk(np.sqrt(Emin[5]), np.insert((-1)*coef[5] , 0, 1))
In [61]:
#Zeros
z
Out[61]:
In [62]:
#poles
p
Out[62]:
In [63]:
angle = np.linspace( 0 , 2 * np.pi , 150 )
radius = 1
x = radius * np.cos( angle )
y = radius * np.sin( angle )
plt.figure(figsize=(10,10))
plt.plot(x, y, 'r', alpha = 0.2)
plt.plot([i.real for i in p], [j.imag for j in p], 'x', label = "Poles")
plt.xlim([-1.1, 1.1])
plt.ylim([-1.1, 1.1])
plt.title("Pole zero plot for LPA of order 6")
plt.xlabel("Real axis")
plt.ylabel("Imaginary axis")
plt.legend()
Out[63]:
Order of lpa p=10¶
In [64]:
z, p, k =tf2zpk(np.sqrt(Emin[9]),np.insert((-1)*coef[9] , 0, 1))
In [65]:
#Zeros
z
Out[65]:
In [66]:
#poles
p
Out[66]:
In [67]:
angle = np.linspace( 0 , 2 * np.pi , 150 )
radius = 1
x = radius * np.cos( angle )
y = radius * np.sin( angle )
plt.figure(figsize=(10,10))
plt.plot(x, y, 'r', alpha = 0.2)
plt.plot([i.real for i in p], [j.imag for j in p], 'x', label = "Poles")
plt.xlim([-1.1, 1.1])
plt.ylim([-1.1, 1.1])
plt.title("Pole zero plot for LPA of order 10")
plt.xlabel("Real axis")
plt.ylabel("Imaginary axis")
plt.legend()
Out[67]:
In [68]:
#can we obtained by sqrt(Emin)
for i in range(0, order):
print(f"Gain for order {i+1} lpa is {np.sqrt(Emin[i])}")
LPC spectra magnitude plots¶
For order p = 2¶
In [69]:
p = 2
signalToEstimate = np.pad(np.array(signalToAnalyse, dtype = float),(p, 3*len(signalToAnalyse)-p), 'constant')
signalEstimate = np.zeros_like(signalToEstimate)
for i in range(0, len(signalToAnalyse) + p-1):
signalEstimate[i] = np.sum(signalToEstimate[i:(i+p)]*coef[p-1])
w, h = freqz(np.sqrt(Emin[p-1]), np.insert((-1*coef[p-1]), 0, 1), None, 1)
plt.figure(figsize=(12, 8)) #plotting the magnitude spectra for the windowed signal
plt.gcf().set_dpi(300)
#plt.subplot(121)
plt.plot(np.arange(0, 8000/2, 8000/dftSize)[0:], 20*np.log10(np.abs(np.fft.fft(signalAfterHamming, dftSize))[:dftSize//2]))
plt.plot(w[:int(len(w)/2)]*8000/(2*3.14), 20*np.log10(np.abs(h[:int(len(h)/2)])))
plt.xlabel("Frequency in Hertz -->")
plt.ylabel("Magnitude of the DFT in dB -->")
plt.legend(["Magnitude of the DFT", f"Magnitude of the LPC of order {p}"])
# plt.subplot(122)
# f = interp1d(find_peaks(20*np.log(np.abs(np.fft.fft(signalToEstimate, dftSize))[:dftSize//2]))[0], 20*np.log(np.abs(np.fft.fft(signalToEstimate, dftSize))[:dftSize//2])[find_peaks(20*np.log(np.abs(np.fft.fft(signalToEstimate, dftSize))[:dftSize//2]))[0]], kind='cubic')
# xNew = np.array([i for i in range(21, 4000)])
# plt.plot(xNew, f(xNew), alpha = 0.3, linewidth = 0.9, color = 'g')
# plt.xlabel("Frequency in Hertz -->")
# plt.ylabel("Magnitude of the DFT in dB -->")
# plt.legend(["Interpolated Magnitude of the DFT"])
Out[69]:
For order p = 4¶
In [70]:
p = 4
signalToEstimate = np.pad(np.array(signalToAnalyse, dtype = float),(p, 3*len(signalToAnalyse)-p), 'constant')
signalEstimate = np.zeros_like(signalToEstimate)
for i in range(0, len(signalToAnalyse) + p-1):
signalEstimate[i] = np.sum(signalToEstimate[i:(i+p)]*coef[p-1])
w, h = freqz(np.sqrt(Emin[p-1]), np.insert((-1*coef[p-1]), 0, 1), None, 1)
plt.figure(figsize=(12, 8)) #plotting the magnitude spectra for the windowed signal
plt.gcf().set_dpi(300)
#plt.subplot(121)
plt.plot(np.arange(0, 8000/2, 8000/dftSize)[0:], 20*np.log10(np.abs(np.fft.fft(signalAfterHamming, dftSize))[:dftSize//2]))
plt.plot(w[:int(len(w)/2)]*8000/(2*3.14), 20*np.log10(np.abs(h[:int(len(h)/2)])))
plt.xlabel("Frequency in Hertz -->")
plt.ylabel("Magnitude of the DFT in dB -->")
plt.legend(["Magnitude of the DFT", f"Magnitude of the LPC of order {p}"])
# plt.subplot(122)
# f = interp1d(find_peaks(20*np.log(np.abs(np.fft.fft(signalToEstimate, dftSize))[:dftSize//2]))[0], 20*np.log(np.abs(np.fft.fft(signalToEstimate, dftSize))[:dftSize//2])[find_peaks(20*np.log(np.abs(np.fft.fft(signalToEstimate, dftSize))[:dftSize//2]))[0]], kind='cubic')
# xNew = np.array([i for i in range(21, 4000)])
# plt.plot(xNew, f(xNew), alpha = 0.3, linewidth = 0.9, color = 'g')
# plt.xlabel("Frequency in Hertz -->")
# plt.ylabel("Magnitude of the DFT in dB -->")
# plt.legend(["Interpolated Magnitude of the DFT"])
Out[70]:
For order p = 6¶
In [71]:
p = 6
signalToEstimate = np.pad(np.array(signalToAnalyse, dtype = float),(p, 3*len(signalToAnalyse)-p), 'constant')
signalEstimate = np.zeros_like(signalToEstimate)
for i in range(0, len(signalToAnalyse) + p-1):
signalEstimate[i] = np.sum(signalToEstimate[i:(i+p)]*coef[p-1])
w, h = freqz(np.sqrt(Emin[p-1]), np.insert((-1*coef[p-1]), 0, 1), None, 1)
plt.figure(figsize=(12, 8)) #plotting the magnitude spectra for the windowed signal
plt.gcf().set_dpi(300)
#plt.subplot(121)
plt.plot(np.arange(0, 8000/2, 8000/dftSize)[0:], 20*np.log10(np.abs(np.fft.fft(signalAfterHamming, dftSize))[:dftSize//2]))
plt.plot(w[:int(len(w)/2)]*8000/(2*3.14), 20*np.log10(np.abs(h[:int(len(h)/2)])))
plt.xlabel("Frequency in Hertz -->")
plt.ylabel("Magnitude of the DFT in dB -->")
plt.legend(["Magnitude of the DFT", f"Magnitude of the LPC of order {p}"])
# plt.subplot(122)
# f = interp1d(find_peaks(20*np.log(np.abs(np.fft.fft(signalToEstimate, dftSize))[:dftSize//2]))[0], 20*np.log(np.abs(np.fft.fft(signalToEstimate, dftSize))[:dftSize//2])[find_peaks(20*np.log(np.abs(np.fft.fft(signalToEstimate, dftSize))[:dftSize//2]))[0]], kind='cubic')
# xNew = np.array([i for i in range(21, 4000)])
# plt.plot(xNew, f(xNew), alpha = 0.3, linewidth = 0.9, color = 'g')
# plt.xlabel("Frequency in Hertz -->")
# plt.ylabel("Magnitude of the DFT in dB -->")
# plt.legend(["Interpolated Magnitude of the DFT"])
Out[71]:
For order p = 8¶
In [72]:
p = 8
signalToEstimate = np.pad(np.array(signalToAnalyse, dtype = float),(p, 3*len(signalToAnalyse)-p), 'constant')
signalEstimate = np.zeros_like(signalToEstimate)
for i in range(0, len(signalToAnalyse) + p-1):
signalEstimate[i] = np.sum(signalToEstimate[i:(i+p)]*coef[p-1])
w, h = freqz(np.sqrt(Emin[p-1]), np.insert((-1*coef[p-1]), 0, 1), None, 1)
plt.figure(figsize=(12, 8)) #plotting the magnitude spectra for the windowed signal
plt.gcf().set_dpi(300)
#plt.subplot(121)
plt.plot(np.arange(0, 8000/2, 8000/dftSize)[0:], 20*np.log10(np.abs(np.fft.fft(signalAfterHamming, dftSize))[:dftSize//2]))
plt.plot(w[:int(len(w)/2)]*8000/(2*3.14), 20*np.log10(np.abs(h[:int(len(h)/2)])))
plt.xlabel("Frequency in Hertz -->")
plt.ylabel("Magnitude of the DFT in dB -->")
plt.legend(["Magnitude of the DFT", f"Magnitude of the LPC of order {p}"])
# plt.subplot(122)
# f = interp1d(find_peaks(20*np.log(np.abs(np.fft.fft(signalToEstimate, dftSize))[:dftSize//2]))[0], 20*np.log(np.abs(np.fft.fft(signalToEstimate, dftSize))[:dftSize//2])[find_peaks(20*np.log(np.abs(np.fft.fft(signalToEstimate, dftSize))[:dftSize//2]))[0]], kind='cubic')
# xNew = np.array([i for i in range(21, 4000)])
# plt.plot(xNew, f(xNew), alpha = 0.3, linewidth = 0.9, color = 'g')
# plt.xlabel("Frequency in Hertz -->")
# plt.ylabel("Magnitude of the DFT in dB -->")
# plt.legend(["Interpolated Magnitude of the DFT"])
Out[72]:
For order p = 10¶
In [73]:
p = 10
signalToEstimate = np.pad(np.array(signalToAnalyse, dtype = float),(p, 3*len(signalToAnalyse)-p), 'constant')
signalEstimate = np.zeros_like(signalToEstimate)
for i in range(0, len(signalToAnalyse) + p-1):
signalEstimate[i] = np.sum(signalToEstimate[i:(i+p)]*coef[p-1])
w, h = freqz(np.sqrt(Emin[p-1]), np.insert((-1*coef[p-1]), 0, 1), None, 1)
plt.figure(figsize=(12, 8)) #plotting the magnitude spectra for the windowed signal
plt.gcf().set_dpi(300)
#plt.subplot(121)
plt.plot(np.arange(0, 8000/2, 8000/dftSize)[0:], 20*np.log10(np.abs(np.fft.fft(signalAfterHamming, dftSize))[:dftSize//2]))
plt.plot(w[:int(len(w)/2)]*8000/(2*3.14), 20*np.log10(np.abs(h[:int(len(h)/2)])))
plt.xlabel("Frequency in Hertz -->")
plt.ylabel("Magnitude of the DFT in dB -->")
plt.legend(["Magnitude of the DFT", f"Magnitude of the LPC of order {p}"])
# plt.subplot(122)
# f = interp1d(find_peaks(20*np.log(np.abs(np.fft.fft(signalToEstimate, dftSize))[:dftSize//2]))[0], 20*np.log(np.abs(np.fft.fft(signalToEstimate, dftSize))[:dftSize//2])[find_peaks(20*np.log(np.abs(np.fft.fft(signalToEstimate, dftSize))[:dftSize//2]))[0]], kind='cubic')
# xNew = np.array([i for i in range(21, 4000)])
# plt.plot(xNew, f(xNew), alpha = 0.3, linewidth = 0.9, color = 'g')
# plt.xlabel("Frequency in Hertz -->")
# plt.ylabel("Magnitude of the DFT in dB -->")
# plt.legend(["Interpolated Magnitude of the DFT"])
Out[73]:
While reconstruction, we would now model the error as white noise since the sound is unvoiced we no longer have the glottal vibrations
In [74]:
signalToEstimate = np.pad(np.array(signalAfterHamming, dtype = float),(0, 10*len(signalAfterHamming)), 'constant')
signalReconstructed = np.zeros_like(signalToEstimate)
mean = 0
std = 1
num_samples = len(signalToEstimate)
whiteNoise = np.random.normal(mean, std, size=num_samples)
plt.figure(figsize = (20, 8))
plt.gcf().set_dpi(300)
plt.plot(whiteNoise)
#plt.vlines(np.arange(0, signalToEstimate.shape[0], 59), 0, 1, colors = 'r', linestyles='--')
plt.xlabel("Times(samples)")
plt.ylabel("Amplitude")
plt.title("White noise generated")
Out[74]:
In [75]:
P = 10#choosing the order of analysis
G = np.sqrt(Emin[P-1])
whiteNoise = np.pad(np.array(whiteNoise, dtype = float),(P, P), 'constant', constant_values=0)
#signal y[n] is assumed casual
for i in range(P, signalToEstimate.shape[0]-2000):
signalReconstructed[i] = G*whiteNoise[i] + np.sum(signalReconstructed[i-P:i]*coef[P-1][::-1])
In [76]:
plt.figure(figsize=(16,8))
plt.gcf().set_dpi(300)
plt.suptitle("Signal estimated along with ground truth")
plt.subplot(121)
plt.plot([i*(1.0/16000) for i in range(len(signalReconstructed))][0:480],signalReconstructed[0:480])
plt.xlabel("Time(seconds)")
plt.ylabel("Amplitude")
plt.title("Estimate of the signal(after resynthesis and de-emphasis)")
plt.subplot(122)
plt.plot([i*(1.0/16000) for i in range(len(signalReconstructed))][0:480],signalToAnalyse[0:480])
plt.xlabel("Time(seconds)")
plt.ylabel("Amplitude")
plt.title("Ground truth(true signal")
Out[76]:
In [77]:
scaled = np.int16(signalReconstructed/np.max(np.abs(signalReconstructed)) * 32767)#scaling the output to 16bit integers for writing to wav file
write('ssReconstructed.wav', 16000, scaled)
In [77]:
In [77]: